Processing seismic data

ABSTRACT

A method of determining a vertical component of particle velocity from seismic data acquired using at least one receiver disposed within a water column comprises determining the vertical component of particle velocity from the pressure at a first location and a second location. The first location is the location of a first receiver, and is vertically beneath the second location. The method uses an operator accurate for seismic data acquired at a vertical separation between the first and second locations of up to at least 0.4 times the minimum wavelength of seismic energy used to acquire the seismic data. The second location may be at or above the surface of the water column. Alternatively, the method be applied to a twin streamer seismic surveying arrangement in which the second location is below the surface of the water column and is the location of a second receiver.

The present invention relates to processing seismic data and, inparticular, to processing marine seismic data to reduce the effects of“ghost reflections”. The present invention can be applied to processingboth marine seismic data acquired in a flat sea and marine seismic dataacquired in a rough sea.

FIG. 1 is a schematic diagram of a marine seismic survey in whichseismic energy is emitted from a source 1 and detected by a seismicreceiver 2 at a depth h below the surface 6 of the sea. Energy emittedfrom the source is reflected by the seabed 3 or by a reflector 4 belowthe seabed 3 and is then detected by the receiver. This path of seismicenergy is labelled 5 in FIG. 1. Information about the geologicalstructure of the earth's interior can be derived from the reflectedseismic energy incident on the receiver The seismic receiver 2 shown inFIG. 1 is a streamer, which is a type of seismic receiver that isfrequently used in marine seismic surveying. A streamer contains aplurality of sensors S₁, S₂, . . . S_(N) such as pressure sensors and/orparticle velocity sensors distributed along its length, which may besome hundreds of metres, and is thus able to measure the reflectedseismic energy at a number of points simultaneously. A streamer may besuspended from one or more floats 8, so that all the receivers of thestreamer are at the same depth in a flat sea.

In addition to the desired path 5 of seismic energy shown in FIG. 1,other seismic energy paths will occur as a result of seismic energybeing reflected or scattered from the sea surface 6. These paths areknown as “ghost reflections”. For example, reference 7 in FIG. 1 shows aghost reflection in which the seismic energy reflected by the reflector4 is not directly incident on the receiver 2, but undergoes a furtherreflection at the sea surface 6 before reaching the receiver. Down-goingsea-surface ghost reflections are an undesirable source of contaminationof seismic data, since they obscure the interpretation of the desiredup-going reflections from the earth's interior.

As is well-known, the ghost signal will produce one or more “notches” inthe frequency spectrum of seismic energy, with the frequency at whichthe notches occur depending on the depth of the receiver below thesea-surface.

The ghost signal reflected from the sea surface is delayed relative tothe direct signal. There are two components to this delay: firstly,there is a 180° phase change upon reflection at the sea surface and,secondly, there is a time delay corresponding to the additional pathlength (which for a signal emitted in the vertical direction is 2z, ortwice the depth of the receiver). The actual vertical far field signalis the sum of the direct signal and the ghost signal. The direct signaland the ghost signal will interfere, and this causes variations in theamplitude of the far field signal. For some frequencies, theinterference is destructive and causes a zero amplitude or “notch” inthe spectrum. This occurs at frequencies where the depth of the receiveris an even number of quarter wavelengths:f _(notch)=(nc/2z), n=0, 1, 2, 3 . . .(where c is the speed of sound in water, n is an integer giving theharmonic order, and z is the depth of the receiver below the seasurface).

Constructive interference occurs at frequencies exactly intermediateadjacent notch frequencies, and this leads to maxima in the amplitude atthese frequencies, given by:f _(peak)=(2n+1)c/4z, n=0, 1, 2, 3 . . .

The effect of the interference between the direct signal and the ghostsignal can be thought of as applying a frequency domain “ghost filter”to the direct signal. The ghost filter has the following form:g(f)=1|+r| ²−2|r|cos(2zf/c)(where r is the reflection coefficient at the sea surface). Theamplitude spectrum of the ghost filter resembles a full-wave rectifiedsine wave, with zeros at the ghost notch frequencies and peaks ofamplitude 2.0 (6 dB) at the peak frequencies.

FIG. 2 shows the amplitude of a typical ghost filter as a function offrequency. This figure shows the ghost filter for the case z=12 m,c=1500 m/s, and with a reflection coefficient of −1.0 at the seasurface. It will be seen that the amplitude decreases to 0 at the notchfrequencies of 0 Hz, 62.5 Hz, 125 Hz . . . , and that there are maximain the amplitude at the peak frequencies of 31.25 Hz, 93.75 Hz . . .

Seismic deghosting is an old and long-standing problem. Traditionally,the primary goal of deghosting has been to broaden the bandwidth of thedata through the notch frequencies. Many prior approaches to de-ghostingseismic data have used the assumptions of a perfectly flat sea surface,a streamer at a known depth below the sea surface, and vanishingpressure at the sea surface. In real seismic data acquisition, however,the sea surface is very often rough. Deterministic algorithms that usean explicit mean streamer depth in the deghosting calculations will notwork reliably in rough sea conditions. Even statistical methods cannotfully correct for the ghost reflections in rough sea conditions, owingto the time variant nature of the problem, as shown by E. Kragh and R.Laws in “Rough seas and statistical deconvolution”, 62^(nd) Annual EAGEMeeting (2001).

The effects of rough seas on the seismic data have recently been thesubject of considerable research. In particular, R. Laws and E. Kraghhave shown, in “Time-lapse seismic and the rough sea wavelet”, 70^(th)Ann. Int. Mtg of Soc. Exploration Geophysicists, Extended Abstracts pp1603-1606 (2000), that errors arising in processing data acquired inrough sea conditions are significant for time-lapse seismic surveyingand for the reliable acquisition of repeatable data for stratigraphicinversion.

As pointed out in UK patent application No 0015810.5, correction of therough-sea surface ghost requires that the wavefield is completelyseparated into its up- and down-going components. The separationrequires in general that the pressure and the vertical component of theparticle velocity (or equivalently the vertical pressure gradient) areboth recorded (see, for example, the equations given by L. Amundsen inGeophysics, Vol 58, pp 1335-1348 (1993)).

Another possibility, as suggested by J. O. A. Robertsson and E. Kragh in“Rough sea deghosting using a single streamer and a pressure gradientapproximation”, submitted to Geophysics (2001) and further investigatedby Rosten et al. in “Rough sea deghosting using a vertical particlevelocity field approximation”, submitted to 63^(rd) Annual EAGE meeting,is to derive an estimate of the vertical component of the particlevelocity from the pressure acquired at a receiver by using the fact thatthe pressure vanishes at the sea surface. The zero pressure along acontour vertically above the streamer may be regarded as a secondstreamer in a dual streamer configuration.

Robertsson and Kragh (above) and PCT application PCT/IB01/01150 suggestusing a technique from synthetic finite-difference modelling of seismicwave propagation, known as the Lax-Wendroff correction, to deriveequations for the vertical pressure gradient at a streamer in thevicinity of a rough sea surface. The equations express the verticalpressure gradient in terms of the pressure gradient along the streamerand the time-derivative of the pressure, both of which may readily beobtained from data acquired by the streamer. The seismic data may thenbe de-ghosted using the pressure and the approximation for the verticalpressure gradient. In contrast to other proposals for seismicdeghosting, the method of Robertsson and Kragh is local: a short spatialfilter (typically of length three) is used to deghost the pressure dataacquired at a single receiver.

A first aspect of the present invention provides a method of determininga vertical component of particle velocity from seismic data acquiredusing at least one receiver disposed within a water column, the methodcomprising: determining the vertical component of particle velocity fromthe pressure at a first location and a second location, the firstlocation being the location of a first receiver and being verticallybeneath the second location, using an operator accurate for seismic dataacquired at a vertical separation between the first and second locationsof up to at least 0.4 times the minimum wavelength of seismic energyused to acquire the seismic data.

The locations may be defined in the co-ordinate system of the firstreceiver, and in this case the locations will be constant with time. Thelocations may alternatively be defined in the earth's co-ordinatesystem, in which case the locations will vary as the first receivermoves through the water.

The present invention provides an improved method of estimating of thevertical particle velocity. Earlier methods of estimating the verticalcomponent of the particle velocity have been applicable only for arestricted range for the depth of the receiver below the surface of thewater column. As discussed below, the prior art method is accurate onlyif the depth of the receiver, below the surface of the water column, isno greater than λ/3.5, where λ is the minimum wavelength of seismicenergy used in the acquisition process. This severely restricts thechoice of receiver arrays with which the method can be used. The presentinvention provides a method of estimating the vertical component of theparticle velocity that can be used with data acquired at greater depthsthan hitherto.

The second location may be at or above the surface of the water column.In this embodiment the invention may be used to process data acquired ata single streamer. In this embodiment, the streeamer may be located at agreater depth within the water column than in the prior art methods.

Alternatively, the second location may be below the surface of the watercolumn and may be the location of a second receiver. In this embodimentthe invention may be used to process data acquired at twovertically-separated streamers. In this embodiment, the verticalseparation between streamers may be greater than in prior art methods.

The method may comprise determining the vertical component of particlevelocity at the second location from pressure acquired at the secondreceiver and from pressure acquired at the first receiver.

The method may comprise determining the vertical component of particlevelocity at the second location according to the following equation orfrom an equation derived therefrom:${v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{{ik}_{z}}{\rho\quad\omega}{\cot\left( {k_{z}\quad\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}} - {\frac{{ik}_{z}}{\rho\quad\omega\quad\sin\quad\left( {k_{z}\Delta\quad z} \right)}{{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}.}}}$

The method may comprise determining the vertical component of particlevelocity at the second location from the following equation:${v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{i}{{\rho\omega\Delta}\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}\quad F_{m}^{(1)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}}} - {\frac{i}{{\rho\omega\Delta}\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(2)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}}}}}$$\begin{matrix}{where} & {{F_{m} = {k^{{- 2}m}\quad{\sum\limits_{n = m}^{\infty}{A_{nm}\left( {k\quad\Delta\quad z} \right)}^{2n}}}},} & {{A_{nm} = {\left( {- 1} \right)^{n - 1}\frac{2^{2n} - 2}{\left( {2n} \right)!}{B_{2n}\begin{pmatrix}n \\m\end{pmatrix}}}},}\end{matrix}$B_(n) is the n^(th) Bernoulli number, and k²=k_(x) ² +k_(y) ².

This provides an approximation to the previous equation.

In a particularly preferred embodiment, each summation over m is carriedout over m=0 to m=1. This provides an approximation that is moreaccurate than known approximations but that may be calculated relativelyeasily.

The operator may be a spatially compact operator. The term “spatiallycompact operator”, as used herein, means a spatial filter having alength (in the spatial domain) that is greater than one but that is lessthan the number of pressure sensors on the streamer.

The method may comprise the further step of processing the seismic datausing the vertical component of the particle velocity. For example, themethod may comprise processing the seismic data so as to determine atleast one of an up-going component and a down-going component of theseismic data.

The method may comprise processing the seismic data using the verticalcomponent of the particle velocity thereby to attenuate effects in theprocessed data of seismic energy reflected and/or scattered at thefree-surface of the water column.

The method may comprise processing the seismic data using the verticalcomponent of the particle velocity thereby to obtain information aboutthe signature of the source of seismic energy.

The seismic data may comprise pressure and particle velocity data, andthe method may comprise comparing the vertical component of the particlevelocity determined from the acquired pressure with the measured valuesof the particle velocity.

Alternatively, the seismic data may comprise pressure and particlevelocity data, and wherein the method may comprise the step ofdetermining the depth of the seismic receiver within the water columnfrom the measured vertical component of the particle velocity and thevertical component of the particle velocity determined from thepressure.

A second aspect of the invention provides a method of determining avertical component of particle velocity from seismic data acquired atfirst and second vertically-separated streamers disposed within a watercolumn, the first streamer comprising M (M>1) pressure sensors disposedalong the length of the first streamer and the second streamercomprising N (N>1) pressure sensors disposed along the length of thesecond streamer, the method comprising applying to the pressure acquiredat sensors on at least one streamer an operator having a length (in thespatial domain) that is greater than one but that is less than thenumber of pressure sensors on the respective streamer.

This embodiment may comprise applying an operator of length m where1<m<M to the pressure acquired at each sensor on the first streamer andapplying an operator of length n where 1<n<N to the pressure acquired ateach sensor on the second streamer.

A third aspect of the invention provides a method of acquiring marineseismic data comprising the steps of: actuating an array of one or moreseismic sources to emit seismic energy; acquiring seismic data at one ormore receivers disposed within a water column, the seismic dataincluding at least pressure data; and processing the acquired pressuredata according to a method as defined above.

A fourth aspect of the invention provides an apparatus for determining avertical component of particle velocity from seismic data acquired usingat least one receiver disposed within a water column, the apparatushaving means for determining the vertical component of particle velocityfrom the pressure at a first location and a second location, the firstlocation being the location of a first receiver and being verticallybeneath the second location, using an operator accurate for seismic dataacquired at a vertical separation between the first and second locationsof up to at least 0.4 times the minimum wavelength of seismic energyused to acquire the seismic data.

A fifth aspect of the invention provides an apparatus for determining avertical component of particle velocity from seismic data acquired atfirst and second vertically-separated streamers disposed within a watercolumn, the first streamer comprising M (M>1) pressure sensors disposedalong the length of the first streamer and the second streamercomprising N (N>1) pressure sensors disposed along the length of thesecond streamer, the apparatus comprising means for determining thevertical component of particle velocity by applying to the pressureacquired at sensors on at least one streamer an operator having a length(in the spatial domain) that is greater than one but that is less thanthe number of pressure sensors on the respective streamer.

The apparatus may comprise a programmable data processor.

A sixth aspect of the invention provides a storage medium containing aprogram for the data processor of an apparatus as defined above.

A seventh aspect of the invention provides a storage medium containing aprogram for controlling a data processor to perform a method as definedabove.

The present invention also provides a method of determining a verticalcomponent of particle velocity from seismic data acquired at a receiverdisposed within a water column, the method comprising determining thevertical component of particle velocity from the pressure acquired atthe receiver using a operator accurate for seismic data acquired at adepth beneath the surface of the water column of up to at least 0.4times the minimum wavelength of seismic energy used to acquire theseismic data. In a preferred embodiment, the operator is accurate forseismic data acquired at substantially any depth beneath the surface ofthe water column.

The operator may be a spatial and/or wavenumber dependent operator.

This method may comprise determining the vertical component of particlevelocity from the pressure acquired at the receiver using the followingequation or an equation derived therefrom: $\begin{matrix}{{v_{z}\left( {\omega,x,y,z} \right)} = {{- \frac{i}{{\rho\omega}\quad z}}{\sum\limits_{m}{{F_{m}\left( {k,z} \right)}\quad{D_{2}^{m}\left( {x,y} \right)}*{p\left( {\omega,x,y,z} \right)}}}}} & (2)\end{matrix}$where p is the acquired pressure wavefield, v_(z) is the verticalcomponent of particle velocity, ω is the angular frequency, k(=ω/α) isthe wavenumber, α is the P-wave velocity, ρ is the density, D₂ is adifferential operator,$F_{m} = {{F_{m}\left( {k,z} \right)} = {k^{{- 2}m}\quad{\sum\limits_{n = m}^{\infty}{A_{nm}({kz})}^{2n}}}}$and * is the 2D spatial convolution operator.

Equation (2) is an exact expression, and is accurate for any receiverdepth. Use of equation (2) enables the vertical particle velocity to beobtained from seismic data acquired at any receiver depth.

In order to reduce the computation required to apply equation (2), it ispossible to solve an equation derived from equation (2) in which thesummation is truncated, rather than solving equation (2) exactly.Solving a truncated form of equation (2) reduces the computationrequired, but at the expense of accuracy. A truncated version ofequation (2) will generally not be accurate for any receiver depth, butwill be accurate only for a receiver depth less than some maximum value.

This method may alternatively comprise determining the verticalcomponent of particle velocity from the pressure acquired at thereceiver using the following equation or an equation derived therefrom:$\begin{matrix}{{{v_{z}\left( {\omega,x,z} \right)} \approx {F_{k}*\quad{p\left( {\omega,x,z} \right)}}},{where}} \\{{\overset{\sim}{F}}_{k} = {\begin{matrix}\min \\{b_{m},a_{n}}\end{matrix}{\sum\limits_{k = {- {({K - 1})}}}^{K - 1}{W_{k}{{\frac{\Theta_{k\quad m}b_{m}}{\Phi_{kn}a_{n}} - F_{k}}}^{2}}}}}\end{matrix}$where a_(n) and b_(m) are the backward and forward filter coefficients,W_(k) are a set of positive weights, * is the 1-D spatial convolutionoperator and F_(k)=F(kΔk_(x)) denote the desired discrete horizontalwavenumber response.

This filter is an approximation of equations (A34) and (A36) given inthe appendix below. Equations (A34) and (A36) are exact, but are notvalid for the case of a rough sea surface; this approximation allowsthem to be applied to the case of a rough sea surface.

This method may use the following filter: $\begin{matrix}{{\overset{\sim}{F}}_{k} = {\begin{matrix}\min \\g_{m}\end{matrix}\quad{\sum\limits_{k = 0}^{\kappa - 1}{W_{k}{{{\Gamma_{km}g_{m}} - F_{k}}}^{2}}}}} \\{{{{where}\quad\Gamma_{km}} = {\cos\left( {k\quad\Delta\quad k_{x}m\quad\Delta\quad x} \right)}},{b_{m} = {0.5g_{m}}}} \\{{{{for}\quad g_{m}} = 1},2,\ldots\quad,{{{M/2}\quad{and}\quad b_{0}} = {g_{0}.}}}\end{matrix}$

The estimate of the vertical particle velocity provided by the presentinvention may be used to de-ghost the data. Alternatively, the estimateof the particle vertical velocity provided by the present invention maybe used in other seismic data processing applications such as estimationof the source signature.

The present invention also provides a method of processing seismic dataacquired at a receiver located on a streamer disposed within a watercolumn, the method comprising the step of determining a derivative ofthe pressure in a horizontal direction perpendicular to the streamerfrom a derivative of the pressure in a direction along the streamer.This may be done using:${\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {\frac{1}{x_{0}}\frac{\partial}{\partial x}{p\left( {x_{0},y_{0}} \right)}}$where the streamer extends along the x-axis and the source is located at(0, y₀).

In this equation, the receiver co-ordinates in the x-y plane are (x₀,y₀). This enables the derivatives of pressure in the cross-linedirection to be estimated even for a single streamer. In contrast, theprior art methods require at least two laterally separated streamers todetermine derivatives in the cross-line direction.

The present invention also provides a method of determining a verticalcomponent of particle velocity comprising determining the verticalcomponent of particle velocity from the pressure acquired at thereceiver using the following equation or an equation derived therefrom:${v_{z}\left( {\omega,x,y,z} \right)} = {{- \frac{i}{{\rho\omega}\quad z}}{\sum\limits_{m}{{F_{m}\left( {k,z} \right)}\quad{D_{2}^{m}\left( {x,y} \right)}*\quad{{p\left( {\omega,x,y,z} \right)}.}}}}$

Other preferred features of the invention are defined in the dependentclaims.

Preferred embodiments of the present invention will now be described byway of illustrative example with reference to the accompanying drawingsin which:

FIG. 1 is a schematic sectional view of a towed marine's seismic survey;

FIG. 2 illustrates the effect of ghost reflections on the amplitudespectrum of seismic energy emitted by the source in the survey of FIG.1;

FIG. 3 is a block flow diagram of a method of the present invention;

FIG. 4 is a schematic plan view of the seismic surveying arrangement ofFIG. 1;

FIGS. 5(a) to 5(d) illustrate the results of a de-ghosting method for astreamer at a depth of 6 m;

FIGS. 6(a) to 6(d) illustrate the results of a de-ghosting method for astreamer at a depth of 10 m;

FIG. 7 is a schematic illustration of a twin streamer seismic surveyingarrangement;

FIGS. 8(a) to 8(d) show results obtained by a twin streamer de-ghostingmethod of the present invention; and

FIG. 9 is a schematic block sectional diagram of an apparatus accordingto the present invention.

Removing the effects of ghost reflections from pressure data acquired ata seismic receiver disposed within a water column is equivalent toseparating the acquired pressure into its up-going and down-goingcomponents. The de-ghosted component of the acquired pressure is theup-going component. Various methods for separating acquired pressuredata into up-going and down-going components have been proposed,including the following: $\begin{matrix}{\overset{\sim}{P} = {\frac{1}{2}\left\lbrack {P - {\frac{\rho\quad\omega}{k_{z}}V_{z}}} \right\rbrack}} & (1)\end{matrix}$

In this equation, {tilde over (P)} denotes the spatial Fourier transformof the de-ghosted (up-going) component of the acquired pressure, Pdenotes the spatial Fourier transform of the acquired pressure p, ρdenotes the density of the water column, ω denotes the angularfrequency, k_(z) is the vertical wave number, and V_(z) is the spatialFourier transform of the vertical component of the particle velocityv_(z).

V_(z), P and {tilde over (P)} will in general be functions of thehorizontal wave numbers k_(x) and k_(y), the angular frequency, and thedepth z of the receiver. The equation (1) above is valid for all wavenumbers.

One prior art method of de-ghosting pressure data acquired at a seismicreceiver disposed within a water column is to use equation (1) todetermine the de-ghosted (up-going) pressure component. This method hashad the disadvantage that, until now, a sufficiently accurate expressionfor v_(z) has not been available. Although various approximateformulations for computing v_(z) have been proposed, these havegenerally only been applicable for a limited range of receiver depth orfrequency of seismic energy. For example, one prior art approximationfor v_(z) is that given in equations (A1) and (A2) in the appendix, andthis prior art approximation is only accurate if the receiver depth isless than approximately λ/3.5, where λ is the minimum wavelength ofseismic energy. This approximation also requires that the separationbetween adjacent receivers on a streamer is around 3 m or less.Furthermore, this prior method is fundamentally limited to frequenciesof seismic energy below the first ghost notch frequency so that, forexample, for a streamer depth of 6 m equation (A1) is accurate up toonly approximately 70 Hz.

The present invention provides more accurate expressions for thevertical particle velocity v_(z), and thus enables more accuratede-ghosting of pressure data using equation (1).

According to the present invention, the vertical component of therequired particle velocity is determined using the following equation,or using an equation derived therefrom: $\begin{matrix}{{v_{z} = {{- \frac{i}{p\quad\varpi\quad z}}\quad{\sum\limits_{m}{{F_{m}\left( {k,z} \right)}\quad{D_{2}^{m}\left( {x,y} \right)}*p}}}},} & (2)\end{matrix}$

In equation (2), the functions F_(m) (k_(z)) are functions that areindependent of the horizontal spatial co-ordinates and are given by:$\begin{matrix}{F_{m} = {{F_{m}\left( {k,z} \right)} = {k^{{- 2}m}\quad{\sum\limits_{n = m}^{\infty}{A_{nm}({kz})}^{2n}}}}} & (3) \\{where} & \quad \\{A_{nm} = {\frac{\left( {- 1} \right)^{n}2^{2n}B_{2n}}{\left( {2n} \right)!}\begin{pmatrix}n \\m\end{pmatrix}}} & (4)\end{matrix}$

In equation (4), A_(nm) are scalar quantities, B_(n) is the n^(th)Bernoulli number, and n! is the factorial of n.

In equation (2) above D₂ is a differential operator, and * denotes thetwo-dimension spatial convolution operator. The components p, v_(z) arein general functions of ω, x, y and z.

Differentiation generally leads to large amplification of highwavenumbers. Hence, in practical situations a bandlimited version of adifferential operator is preferably used for D₂.

The expression for v_(z) given in equation (2) is valid for all receiverdepths in the water column. Use of equation (2) to determine thevertical component of the particle velocity enables accuratedetermination of the vertical component of the particle velocity forpressure data acquired at any depth in a water column. This in turnenables further processing steps that require knowledge of the verticalcomponent of the particle velocity—such as de-ghosting, for example—tobe carried out accurately for pressure data acquired at any depth in awater column.

The present invention is not limited to use of the exact equation (2) toobtain the vertical component of the vertical particle velocity. Theinvention also envisages that the vertical component of the verticalparticle velocity may be obtained using an equation derived fromequation (2) such as, for example, an approximate solution to equation(2).

A special case of equation (2) arises when m is finite, and extends fromzero to M. In this case, and assuming infinite bandwidth, equation (2)may be re-written as: $\begin{matrix}{{v_{z} \approx {{- \frac{\mathbb{i}}{p\quad\varpi\quad z}}{\sum\limits_{m = 0}^{M}{{F_{m}\left( {k,z} \right)}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)^{m}\quad p}}}},} & (6)\end{matrix}$

The functions F_(m) are given in equation (3) above as infinite series.The functions F_(m) may also be written analytically in terms oftrigonometric functions, and examples for F₀, F₁, F₂ are given inequations (A8), (A9) and (A10) below.

It will be seen that the x- and y-variables in equation (6) arede-coupled from z. This indicates that equation (6) is valid for a roughsurface of the water column, as well as for a flat surface of the watercolumn. This preferred embodiment of the present invention thereforeprovides a further advantage, since it enables accurate de-ghosting ofpressure data acquired in rough sea conditions to be carried out, forany receiver depth in the water column.

In the case where M=2 equation (6) reduces to: $\begin{matrix}{v_{z} \approx {{- {\frac{\mathbb{i}}{\rho\quad\omega\quad z}\left\lbrack {F_{0} + {F_{1}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)} + {F_{2}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)}^{2}} \right\rbrack}}p}} & (7)\end{matrix}$

The spatial derivatives in equations (6) and (7) may be estimated fromthe pressure acquired at two or more receivers. For example, thepressure derivative along the x-axis (which is assumed to coincide withthe axis of the streamer) may be determined from the pressure recordedat the same time at different receivers along the streamer.

The preferred method of determining the spatial derivatives will dependon the spatial sampling interval—that is, the distance between adjacentsensors on a streamer. Where there is a low sampling interval—that isthe sensors are closely spaced and so provide densely sampled pressuredata, the spatial derivatives may be represented by classical centralfinite—difference approximations. In the case of equation (7), forexample, the second-order spatial derivatives may be estimated using athree-point finite difference operator that estimates the spatialderivative in the x-direction of the pressure at one receiver from thepressure values acquired at that receiver and at the adjacent receiveron each side (this is sometimes known as an iterative filter of length3).

For a coarse receiver spacing, the horizontal derivatives are mostpreferably implemented as numerically optimised difference operators, orusing Fourier techniques. In a Fourier technique, the derivatives areobtained by Fourier transforming the pressure to the wavenumber domain,multiplying by the negative square of the wavenumber (for the secondderivative), and inverse Fourier transforming the result.

Equations (6) and (7) are derived from, and are particular cases of,equation (2), in which the variable m is chosen to be finite. Furtherspecial cases occur when the variables m and n in equation (2) are bothchosen to be finite, with m varying from 0 to M and with n varying fromm to N (m). In this case, equation (2) reduces to the following:$\begin{matrix}{v_{z} \approx {\frac{\mathbb{i}}{p\quad\varpi\quad z}{\sum\limits_{m = 0}^{M}{{F_{m}\left( {k,z} \right)}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)^{m}\quad p}}}} & (8)\end{matrix}$

In equation (8), the functions F_(m) have the same general form as givenin equation (3), but the summation over n is now carried out from m toN(m), as shown in equation (A20) below.

If we choose M=1, N(0)=2, and N(1)=1, equation (8) reduces to thefollowing: $\begin{matrix}{v_{z} \approx {{- {\frac{\mathbb{i}}{p\quad\varpi\quad z}\left\lbrack {A_{00} + {A_{10}({kz})}^{2} + {A_{11}{z^{2}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)}}} \right\rbrack}}\quad p}} & (9)\end{matrix}$where the scalar coefficients have the following values: A₀₀=1,A₁₀=A₁₁=−⅓.

It will be seen that equation (9) has the same general structure as theprior art equation (A1) given in the Appendix below. The scalarcoefficients A₁₀ and A₁₁ in equation (9), however differ from the scalarcoefficients in equation (A1). This has the significant effect thatequation (9) is valid for greater streamer depths than equation (A1).(Equation (2) above is valid for all streamer depths, but the depthrange over which an approximate solution to equation (2) will be validvaries from one approximate solution to another.) Equation (9) is validat depths up to approximately λ/2.5, where λ is the minimum wavelength,whereas equation (A1) requires that the streamer depth is no greaterthan approximately λ/3.5.

FIG. 3 is a block flow diagram illustrating one method of the presentinvention.

Initially, at step 9, seismic data are acquired at a plurality ofreceivers disposed within a water column. In general, step 9 willconsist of acquiring seismic data using a seismic surveying arrangementof the general type shown in FIG. 1 in which a plurality of seismicreceivers S1 . . . SN are disposed on a streamer 2 that is set at adepth h below the surface of the water column. The data may be acquiredwith a single streamer, or they may be acquired with two or morestreamers that are separated laterally (that is, separated along they-axis) from one another.

The invention may alternatively be used to process pre-existing seismicdata Step 9 of acquiring seismic data may therefore be replaced by step10 of retrieving pre-existing seismic data from storage.

At step 11, the vertical component of the particle velocity v_(z) iscalculated for one receiver (for example the j^(th) receiver). Step 11is carried out using equation (2), or using any equation derived fromequation (2) such as equation (6), (7), (8) or (9). As outlined above,the values of the pressure derivatives occurring in these equations areestimated using pressure values acquired at two or more receivers.

Once the vertical component of the particle velocity has been determinedfor the j^(th) receiver, the pressure data acquired at the j^(th)receiver may then be de-ghosted at step 12. The de-ghosting step iscarried out using equation (1), with the component V_(z) being obtainedby a spatial Fourier transform on the output of step 11.

At step 13 the results of step 12 are output, for example by display byan operator, or are alternatively stored for subsequent retrieval.

A receiver S_(i) on the streamer 2 of FIG. 1 will acquire data bysampling the incident seismic wavefield at discrete time values and sowill record the values of the wavefield at times t₀, t₀+t_(s), t₀+2t_(s) . . . , where t_(s) is the sampling time. Steps 11 to 13 have beenperformed on data acquired at one particular sampling time. At step 14,therefore it is determined whether the de-ghosting should be repeatedfor data acquired by the receiver at a different sampling time. If thisyields a “yes” determination, a new sampling time is selected at step15, and steps 11-13 are repeated for data acquired at the new time. Thisprocess is repeated until a “no” determination is obtained at step 14.

At step 16 it is determined whether steps 11, 12 and 13 should berepeated for data acquired by another receiver on the streamer. If step16 yields a “yes” determination, a new receiver is selected at step 17,and steps 11 to 14 and, if necessary step 15, are repeated for the newreceiver. This loop is repeated until a “no” determination is obtainedat step 16.

If it is desired systematically to de-ghost data acquired by everyreceiver on a streamer, the index j may initially be set to 1, step 16may consist in determining whether the index j is equal to the totalnumber of receivers J on the receiver, and step 17 may consist ofincrementing j by 1.

When a “no” determination is obtained at step 16, step 18 thendetermines whether the process should be repeated for data acquired byanother streamer. If this step yields a “yes” determination, a newstreamer is selected at step 19, and steps 11-17 are then repeated forthe receivers on the new selected streamer.

Where the method of the present invention is applied to data acquiredusing only a single streamer, step 18 may be omitted and step 19 is notapplicable.

As noted above, the present invention may be applied to seismic dataacquired using only a single streamer—that is, to seismic data acquiredat receivers all located at the same y-co-ordinate value. In this case,it would initially appear difficult to obtain the spatial derivatives inthe y-direction (perpendicular to the length of the streamer), that arerequired to allow solution of, for example, equations (6), (7), (8) or(9). A further aspect of the present invention, therefore, provides amethod of estimating the necessary derivatives in the y-direction aderivative of the pressure along the streamer. A streamer generallycontains a large number of pressure sensors disposed at regularintervals along the length of the streamer, and the spatial derivativeof the pressure along the streamer can readily be estimated from thepressure acquired at the sensors. The present invention makes itpossible to estimate a spatial derivative in a direction perpendicularto the streamer from data acquired for a single streamer, from a spatialderivative of the pressure along the streamer.

The method of estimating a spatial derivative in a directionperpendicular to the streamer from data acquired for a single streamermay be carried out in connection with a determination of the verticalcomponent of the particle velocity according to the method describedabove, but it is not limited to use in connection with determination ofthe vertical component of the particle velocity.

In a preferred embodiment, a spatial derivative in a directionperpendicular to the streamer may be estimated using the followingequation: $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {\frac{1}{x_{0}}\frac{\partial}{\partial x}{{p\left( {x_{0},y_{0}} \right)}.}}} & (10)\end{matrix}$

The derivation of equation 10 is explained in detail in the Appendixbelow.

FIGS. 5(a) to 5(d) and 6(a) to 6(d) illustrates the results provided bythe present invention. These results were obtained by modellingsynthetic seismic data for a rough-sea surface, and de-ghosting thesynthetic seismic data using the method of the present invention. Thesynthetic seismic data were computed using the rough-sea modellingmethod developed by R. Laws and E. Kragh (above) and by E. Kragh and R.Laws (above). In summary, rough-sea surfaces are generated, and impulseresponses beneath these surfaces are computed. The impulses are made insets of three, and a total of 40 sets were computed with each set beingcomputed from a different realisation of a 4 m significant wave heightsea surface. Data were synthesised for a streamer depth of 6.0 m (FIGS.5(a) to 5(d)) and a streamer depth of 10.0 m (FIGS. 6(a) to 6(d)).

The synthetic seismic data were then de-ghosted using equation (7) aboveto compute the vertical particle velocity. In the computation usingequation (7), the parameters were set as follows: M=1, N(0)=∞ andN(1)=2, so giving F₀=kz cot (kz), $\begin{matrix}{{F_{1} = {{- \frac{z^{2}}{3}}\left( {1 + {\frac{2}{15}k^{2}z^{2}}} \right)}},} & {{{and}\quad F_{2}} = 0.}\end{matrix}$

A zero-phase Ricker wavelet with a maximum frequency of 50 Hz wasconvolved with the 120 pulses for estimating the v_(z) component. In theestimation of the v_(z) component, the well-known length-three centraldifference operator was used to obtain the second-order horizontalspatial derivative in the x-direction. (The calculation was performedusing a 2-D estimation, so that derivatives in the cross-line direction(y-direction) become zero and may be ignored.)

Once the v_(z) component had been estimated using equation (7), thede-ghosted pressure was determined using equation (1) above. Themodelled pressure and the estimated v_(z) component were down-sampled bya factor of 4 in time before the de-ghosted pressure was determined,thereby increasing the temporal sampling distance from 0.25 ms to 1.0ms.

FIG. 5(a) shows the raw synthetic seismic data, indicated generally asA, and the results of de-ghosting the synthetic data using theLax-Wendroff technique of Robertsson and Kragh (above), indicated as B.

FIG. 5(b) shows the raw synthetic seismic data, indicated generally asA, and the results of de-ghosting the synthetic data using equation (7)indicated as C. FIG. 5(c) shows the 95% confidence interval in spread ofthe amplitude spectrum, and FIG. 5(d) shows the 95% confidence intervalin the spread of the phase spectrum. The results are shown for the rawsynthetic data, indicated as A, the results of de-ghosting the syntheticdata using the prior art Lax-Wendroff technique, indicated as B, and theresults of de-ghosting the synthetic data using equation (7) indicatedas C. Although the rough-sea perturbations are not completely removedfrom the de-ghosted data obtained by the present invention, thecorrection shows a clear improvement below the first ghost-notchfrequency (approximately 125 Hz). The residual error is of the order of0.5-1.0 dB in amplitude and about 5° in phase, and these errors aresignificantly lower than the errors currently obtained in de-ghostingreal seismic data.

FIGS. 6(a) to 6(d) correspond generally to FIGS. 5(a) to 5(d)respectively, except that they relate to a streamer depth of 10.0 m. Fora streamer at this depth, the first ghost notch frequency isapproximately 75 Hz.

A zero-phase band pass filter with cut off frequencies of 10-15-90-100Hz was applied to the wavelets shown in FIG. 2, whereas a zero-phaseband pass filter with cut-off frequencies 10-15-65-70 was applied to thewavelets shown in FIG. 6(a) to 6(d). This has a similar effect to thebandlimited differentiation mentioned above and also below in connectionwith equations (A14) and (A15). Here the filter is applied to the datadirectly in a pre-processing step, and this has a similar overalleffect.

The de-ghosted seismic data obtained by a method of the presentinvention may be subjected to further processing. For example, thede-ghosted data may be used to estimate the signature of the seismicsource 1. This can be done, for example, according to the methodproposed by L. Amundsen et al in Geophysics Vol 60, pp 212-222 (1995).

Alternatively, the de-ghosted data may be further processed to remove orat least attenuate events in the data arising from paths that involvereflection and/or scattering at the free-surface of the water column.Such events are generally known as free-surface related multiples, sincethey arise from paths that make multiple passes through the watercolumn. In one embodiment, the principal steps required to removefree-surface related multiples are as follows.

Initially the direct arrival and the ghost arrival associated with thedirect arrival in the incident pressure wavefield acquired at thestreamer are isolated and removed from the pressure wavefield. This maybe done in any suitable way, for example, by “mute” in which the directarrival and its associated ghost arrival is simply separated from theother arrivals. Provided that the source-receiver offset is not toogreat and the water depth is not too shallow, the direct arrival in theseismic data, and its associated ghost arrival, are generally wellseparated from other arrivals.

The vertical component of the particle velocity is then estimated,according to a method as described above. The vertical component of theparticle velocity is then used to separate the remaining incidentpressure wavefield (i.e., the pressure wavefield remaining after thedirect arrival and its associated ghost have been removed) into itsup-going and downgoing components, for example using equation (1) aboveto determine the up-going component of the pressure wavefield.

Once the downgoing component of the remaining part of the pressurewavefield has been found, it is added to the original direct arrival andits associated ghost arrival. This gives the overall downgoingwavefield. (In a typical survey the source will be nearer the seasurface than the receiver, so that the direct arrival and its associatedghost contain only down-going energy.) Effects arising from free-surfacemultiple reflections may then be removed using the “U divided by D”method proposed by L. Amundsen in “Elimination of free-surface relatedmultiples without need of the source wavelet”, Geophysics, Vol. 66, pp327-341 (2001). In this method, “U divided by D” refers to the up-goingcomponent of the wavefield remaining after removal of the direct arrivaland its associated ghost arrival, divided by the overall down-goingcomponent of the wavefield. This division is essentially a deconvolutionstep.

The embodiments of the invention described above may be carried out onseismic data acquired by one or more pressure sensors disposed within awater column, for example on pressure data acquired by a streamer.Alternative embodiments of the invention will now be described in whichthe seismic data include both pressure data and particle velocity data.Such data may be acquired by, a seismic receiver, for example astreamer, that has both pressure sensors and particle velocity sensors.A streamer provided with both pressure sensors and particle velocitysensors can directly measure both the fluid pressure p and the particlevelocity v.

The quantities p and v are not independent from one another, but arerelated by the following equations: $\begin{matrix}{\overset{.}{p} = {{- K}{\nabla{\cdot v}}}} & (11) \\{v = {{- \frac{1}{\rho}}{\nabla p}}} & (12)\end{matrix}$

In equation (11), {dot over (p)} denotes the time-derivative of p, and Kdenotes the bulk modulus of the water column.

In one embodiment of the invention, it is assumed that the depth of thestreamer is known. The depth may be known if, for example, the streameris provided with depth sensors. If the streamer depth is known, it isthen possible to compare the output from the particle velocity sensorswith the output from the fluid pressure sensors. This comparison can bemade, for example, by calculating the particle velocity from themeasured values for pressure and from the depth of the streameraccording to equation (2) or an equation derived therefrom, andcomparing the calculated value of the particle velocity with directlymeasured values for the particle velocity. These values should, ofcourse, be the same if the particle velocity and the fluid pressure havebeen measured accurately. This embodiment can be considered ascalibrating the pressure sensors relative to the particle velocitysensors, and allows the accuracy of the sensors to be monitored.

This process for monitoring the accuracy of the sensors can be carriedout while acquired seismic data is being processed by a method of thepresent invention. Alternatively, it can be carried out on its own, forexample to check the sensors at the start of a seismic survey. Thisallows the extra degree of freedom provided by measuring both p and v tobe used to calibrate the sensors.

In this embodiment, the depth of the streamer can be estimated, forexample from knowledge of the arrangement for suspending the streamerfrom floats. For improved accuracy, however, it is preferable that thedepth of the receiver is measured directly.

The invention may also be used to determine the depth of the receiver.To do this the vertical component of the particle velocity is determinedusing a method as described above. The depth of the receiver may beestimated from the measured value for the vertical component of theparticle velocity and the estimated value for the vertical component ofthe particle velocity.

In an alternative embodiment of the invention, the vertical component ofparticle velocity from the pressure acquired at the receiver isdetermined using the following equation or an equation derivedtherefrom: $\begin{matrix}{{{v_{z}\left( {\omega,x,z} \right)} \approx {F_{k}*{p\left( {\omega,x,z} \right)}}},} & (13) \\{{{where}\quad{\overset{\sim}{F}}_{k}} = {\min\limits_{b_{m},a_{n}}{\sum\limits_{k = {- {({K - 1})}}}^{K - 1}{W_{k}{{\frac{\Theta_{km}b_{m}}{\Phi_{kn}a_{n}} - F_{k}}}^{2}}}}} & (14)\end{matrix}$

In equation (14), a_(n) and b_(m) denote the backward and forward filtercoefficients, W_(k) are a set of positive weights, and F_(k)=F(kΔk_(x))denote the desired discrete horizontal wavenumber response. Thederivation of equation (14) is explained in the Appendix below.

As noted above, equation (14) is an approximation of equations (A34) and(A36) given in the appendix below. Solutions to equation (14) exist forany depths. Equation (14) may be applied to data obtained with a roughsea surface, which is not possible for equations (A34) and (A36)(although they are exact for the case of a flat sea surface). Asexplained in the Appendix, equation (14) may be implemented to give thefollowing zero-phase non-recursive filter: $\begin{matrix}\begin{matrix}{{\overset{\sim}{F}}_{k} = {\min\limits_{g_{m}}\quad{\sum\limits_{k = 0}^{\kappa - 1}{W_{k}{{{\Gamma_{km}g_{m}} - F_{k}}}^{2}}}}} \\{{{{where}\quad\Gamma_{km}} = {\cos\left( {k\quad\Delta\quad k_{x}m\quad\Delta\quad x} \right)}},{b_{m} = {0.5g_{m}}}} \\{{{{for}\quad g_{m}} = 1},2,\ldots\quad,{{{M/2}\quad{and}\quad b_{0}} = {g_{0}.}}}\end{matrix} & (15)\end{matrix}$

Vertical particle velocities obtained using equations (13) and (14) orusing equations (13) and (15) may be used in any of the applicationsdescribed herein. Equations (13), (14) and (15) are 1-D equations, butthe corresponding 2-D equations may also be used. For the 2-D case, theF_(k) are a function of both k_(x) and k_(y).

The embodiments of the present invention described above have related tothe de-ghosting of the signal recorded at the receiver. The inventioncan, however, also be used for source-side de-ghosting. In principlethis could be done by directly measuring the fluid pressure at thesource, but in practice it is possible to make use of the principle ofreciprocity to avoid the need to do this.

The principle of reciprocity is a fundamental principle of wavepropagation, and states that a signal is unaffected by interchanging thelocation and character of the sources and receivers. For example, if asurveying arrangement with a seismic source at point A and a receiverarray at point B gives a certain signal at the receiver array, thenusing a single receiver at point A and a source array at point B wouldlead to the same signal, provided that the source array corresponds tothe receiver array. (By “corresponds”, it is meant that the source arraycontains the same number of sources as the receiver array has receivers,and that the sources in the source array are arranged in the samepositions relative to one another as the receivers in the receiverarray).

In consequence of the principle of reciprocity, the signal emitted by anarray of seismic sources can be de-ghosted by making measurements usinga receiver array, whether the seismic sources are actuatedsimultaneously or consecutively. The recorded signal is analogous tothat recorded in a reciprocal arrangement where the source and receiverlocations are interchanged. The method outlined hereinabove cantherefore also be applied to the source array that it is desired tode-ghost. The signal produced at the receiver array by an array of aplurality of seismic sources is measured, and a de-ghosting filter isderived from this measured signal. By the principle of reciprocity, thisfilter can be used to de-ghost the signal emitted by the source array.

The methods described in this application are fast, and it is possibleto process a complete trace using a method of the invention. If desired,however, the methods may be applied to a selected part of a trace.

FIG. 7 shows a twin streamer seismic surveying arrangement in which twostreamers are disposed with a water column at different depths below thesurface of the water column. Components having the same references inFIG. 7 as in FIG. 1 correspond to the components described above withreference to FIG. 1, and the description of these components will not berepeated.

In the seismic surveying arrangement of FIG. 7, a first streamer 2 issuspended within the water column. A second streamer 2′ is suspendedbelow the first streamer 2. The second streamer 2′ is deployed so as tobe nominally vertically below the first streamer 2, although the actionof tides and currents may mean that the second streamer is not alwaysexactly vertically below the first streamer. The streamers are deployedsuch that the first streamer is deployed at a nominal depth z, below thesurface 6 of the water column and the second streamer is deployed at anominal depth z₂ (where Z₂>Z₁) below the surface 6 of the water column.Again, the action of tides and currents, and the effect of surfacewaves, may mean that the first and second streamers are not always atexactly their nominal depths; the separation between the streamers mayalso vary form the intended value of Δz (where Δz=z₂−Z₁).

A plurality of pressure sensors is disposed along each streamer, withthe upper streamer 2 having pressure sensors S1, S2 . . . SN up to atotal of N pressure sensors and the lower streamer 2′ having pressuresensors S′1, S′2 . . . S′M up to a total of M pressure sensors. Eachstreamer may have the same number of pressure sensors (in which caseN=M). One or both of streamers may also be provided with other sensorssuch as, for example, particle velocity sensors.

The seismic surveying arrangement also comprises an apparatus 20 forprocessing seismic data acquired by the sensors on the streamers 2,2′ asa result of actuating the source array 1 to emit seismic energy. Theprocessing apparatus may be located, for example, on shore, on thetowing vessel, or on another vessel. Data acquired at the sensors on thestreamers may for example be transmitted to a storage means (for examplelocated on the towing vessel) by, for example, an electrical, optical orwireless link (not shown), and may subsequently be passed to theprocessing apparatus. Alternatively, data acquired at the sensors may bestored for subsequent retrieval in storage means provided on thestreamers.

There is a known method of estimating the up-going pressure at the upperstreamer from the pressure measurements acquired at the pressure sensorson the two streamers (see equation A41 of the appendix). However, thisprior method is limited to the case of horizontal streamers having aconstant vertical separation.

The method of the present invention makes it possible to obtain anestimate of the vertical component of the particle velocity at the upperstreamer from the pressure measurements acquired at the pressure sensorson the two streamers. The vertical component of the particle velocity atthe upper streamer may be found from the acquired pressure measurementsaccording to: $\begin{matrix}{{v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{{ik}_{z}}{\rho\quad\omega}{\cot\left( {k_{z}\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}} - {\frac{{ik}_{z}}{\rho\quad\omega\quad\sin\quad\left( {k_{z}\Delta\quad z} \right)}{{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}.}}}} & (16)\end{matrix}$

As explained in the Appendix, equation (16) may be approximated as:$\begin{matrix}{{v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{i}{\rho\quad\omega\quad\Delta\quad z}{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(1)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}}} - {\frac{i}{{\rho\omega\Delta}\quad z}{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(2)}\kappa^{2m}{{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}.}}}}}} & (17)\end{matrix}$

In the simplest approximation of equation (17), only the lowest term isretained in each summation—that is, F₀ ⁽¹⁾ in the first summation and F₀⁽²⁾ in the second summation. This is equivalent to the known “shift andsum” technique. According to the invention, however, the particlevelocity is estimated in a higher order approximation, in which at leastone more term is retained in at least one of the summations. It ispossible to obtain an estimate of the vertical particle velocity at thelower streamer or at other locations that do not coincide with eitherstreamer in a similar way using equivalent equations readily apparent tothose skilled in the art.

In a particularly preferred-embodiment of the invention, the first twoterms are retained in each summation, so that F₀ ⁽¹⁾ and F₁ ⁽¹⁾ areretained in the first summation and F₀ ⁽²⁾ and F₁ ⁽²⁾ are retained inthe second summation. This can be implemented with a 3-point runningfilter.

FIGS. 8(a) to 8(d) show results for the 3-point filter approximation inwhich the first two terms are retained in each summation of equation(17). FIGS. 8(a) to 8(d) show how the simple “shift & sum” approximationand the more accurate 3-point filter approximation compare to the truesolution for estimating vertical particle velocity in a twin-streamerconfiguration as given by equation (16). Since equations (16) and (17)contain two terms with one filter being applied to the shallowerstreamer data and one filter being applied to the deeper streamer data,FIGS. 8(a) to 8(d) show the results for the two terms separately.

FIGS. 8(b) and 8(d) show the difference between the “true” solution forestimating the vertical component of the particle velocity from equation(16) (as obtained using a numerical technique) and the “shift & sum”approximation of equation (17) using just the first term in eachsummation. The results were obtained for a constant streamer separationof 7 m. FIG. 8(b) shows the difference between the “true” solution andthe “shift and sum” approximation for the first term (the term in thepressure at the upper streamer) and FIG. 8(d) shows the differencebetween the “true” solution and the “shift and sum” approximation forthe second term (the term in the pressure at the lower streamer).

FIGS. 8(a) and 8(c) show the difference between the “true” solution forestimating the vertical component of the particle velocity from equation(16) using a numerical technique and the 3-point filter approximation ofequation (17) using the first two terms in each summation. The resultswere again obtained for a constant streamer separation of 7 m. FIG. 8(a)shows the difference between the “true” solution and the 3-point filterapproximation for the first term (the term in the pressure at the upperstreamer) and FIG. 8(c) shows the difference between the “true” solutionand the 3-point filter approximation for the second term (the term inthe pressure at the lower streamer).

The differences shown in FIGS. 8(a) and (b) have been normalised againstthe absolute value of the first term (the term in the pressure at theupper streamer) in the true solution. The differences shown in FIGS.8(c) and (d) have been normalised against the absolute value of thesecond term (the term in the pressure at the lower streamer) in the truesolution.

It will be seen that, although the “shift & sum” algorithm is areasonable approximation, especially at low wave numbers, the 3-pointfilter approximation provides better results.

The regions of FIGS. 8(a) and 8(c) in which the difference between the3-point filter approximation and the true solution is significant areindicated in the Figures and are considerably smaller than the regionsof FIGS. 8(b) and 8(b) in which the difference between the shift and sumapproximation and the true solution is significant.

A further advantage of the invention is that it may be applied even ifthe vertical separation between the streamers varies over the length ofthe streamers. It is not limited to the case of a constant separationbetween the streamers.

As explained above, if the upper streamer is provided with particlevelocity sensors it is possible to compare the measured value for v_(z)with the value for v_(z) obtained from equation (17), for example byusing the 3-point filter approximation. This allows the depth of thestreamer to be determined, or allows the particle velocity sensors to becalibrated as explained above.

As also explained above, the value for v_(z) obtained from equation(17), for example by using the 3-point filter approximation, may be usedin de-ghosting the seismic data acquired at the streamers or it may beused to eliminate effects of the source signature.

A further application of the twin streamer seismic surveying arrangementof FIG. 7 is in the use of a Lax-Wendroff technique to obtain thevertical component of the particle velocity from pressure acquired atthe streamers. The vertical component of the particle velocity isequivalent to the vertical pressure gradient ∂p/∂z, and the Lax-Wendrofftechnique allows the vertical pressure gradient ∂p/∂z to be obtainedfrom the time derivative of the pressure and from ∂p/∂x which is thederivative of the pressure along the streamer. The Lax-Wendrofftechnique is described in PCT/IB01/01150, the contents of which areincorporated by reference.

PCT/IB01/01150 gives the following equation for the vertical pressuregradient at a point (x₀,z₀) below the surface of a water column:$\begin{matrix}{{\partial_{z}{p\left( {x_{0},z_{0}} \right)}} = {\frac{{p\left( {x_{0},z_{0}} \right)} - {p\left( {x_{0},{z_{0} - h}} \right)}}{h} + {\frac{h}{2}{\partial_{zz}{p\left( {x_{0},z_{0}} \right)}}} + {O\left( h^{2} \right)}}} & (18)\end{matrix}$

In equation (18), h denotes the vertical distance of the streamer at apoint (x₀,z₀) below the surface of the water column. PCT/IB01/01150simplifies equation (18) by using the fact that the pressure vanishes atthe surface of the water column so that p(x₀,z₀−h) is zero, and byreplacing the term in ∂_(zz)p(x₀,z₀) according to: $\begin{matrix}{{\partial_{zz}{p\left( {x_{0},z_{0}} \right)}} = {{\frac{1}{\alpha^{2}}{\overset{¨}{p}\left( {x_{0},z_{0}} \right)}} - {\partial_{zz}{p\left( {x_{0},z_{0}} \right)}}}} & (19)\end{matrix}$

In the twin streamer case, pressure data is available at two depths.Equation (18) may be applied to the twin streamer case, by setting hequal to Δz(=z₁−z₂). The pressure acquired at the pressure sensors onthe two streamers provides the two terms in p in equation (18), and theterm in ∂_(zz)p may be estimated from the pressure acquired at differentsensors along a streamer.

Thus, the invention provides a further method of determining thevertical component of particle velocity from pressure measurements madeusing a twin streamer arrangement. In this method pressure measurementsacquired at n pressure sensors on the upper streamer 2, where 1<n<N, andpressure measurements acquired at m pressure sensors on the lowerstreamer 2′, where 1<m<M, are used in the determination of the verticalcomponent of particle velocity.

FIG. 9 is a schematic block diagram of an apparatus 20 that is able toperform a method according to the present invention.

The apparatus 20 comprises a programmable data processor 21 with aprogram memory 22, for instance in the form of a read only memory (ROM),storing a program for controlling the data processor 21 to processseismic data by a method of the invention. The apparatus furthercomprises non-volatile read/write memory 23 for storing, for example,any data which must be retained in the absence of a power supply. A“working” or “scratch pad” memory for the data processor is provided bya random access memory RAM 24. An input device 25 is provided, forinstance for receiving user commands and data. One or more outputdevices 26 are provided, for instance, for displaying informationrelating to the progress and result of the processing. The outputdevice(s) may be, for example, a printer, a visual display unit, or anoutput memory.

The seismic model, the parameters of a seismic surveying arrangement andthe prior AVO information may be supplied via the input device 25 or mayoptionally be provided by a machine-readable data store 27.

The results of the processing may be output via the output device 26 ormay be stored.

The program for operating the system and for performing the methoddescribed hereinbefore is stored in the program memory 22, which may beembodied as a semiconductor memory, for instance of the well known ROMtype. However, the program may well be stored in any other suitablestorage medium, such as a magnetic data carrier 22 a (such as a “floppydisk”) or a CD-ROM.

Appendix

Robertsson and Kragh (above) have shown that the vertical component ofthe particle velocity, v_(z), can be estimated from the raw measurementsof the pressure field, p. In the typical seismic frequency band onepossible approximation given is $\begin{matrix}{{{v_{z}\left( {\omega,x,y,z} \right)} = {{- {\frac{i}{\rho\quad\omega\quad z}\left\lbrack {A_{00}^{\prime} + {A_{10}^{\prime}({kz})}^{2} + {A_{11}^{\prime}{z^{2}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)}}} \right\rbrack}}\quad{p\left( {\omega,x,y,z} \right)}}},} & ({A1}) \\\begin{matrix}{{with}\quad{scalars}} & {{A_{00}^{\prime} = 1};} & {{A_{10}^{\prime} = {{- 1}/2}};} & {{A_{11}^{\prime} = {{- 1}/2}},}\end{matrix} & ({A2})\end{matrix}$

In equation (A1) ω is the (angular) frequency, k=ω/α is the wavenumber,α is the P-wave velocity, ρ is the density, and z=z(x,y) is the streamerdepth. According to Robertsson and Kragh, approximation requires thatthe streamer is towed no deeper than approximately λ/3.5, where λ is theminimum wavelength, and that the receiver spatial sampling interval isaround 3 m or less. This method is fundamentally limited to frequenciesbelow the first ghost notch. For instance, for a streamer depth of 6 m,equation (A1) is accurate up to approximately 70 Hz.

For densely sampled pressure measurements one attractive feature ofmethod (A1) is that the filters for estimating the vertical particlevelocity are local and short. By approximating the second-orderhorizontal spatial derivatives by three-point centraldifference-operators, the spatial filter has cross structure with fivecoefficients. Processing single streamer pressure data such as thoseprovided by the new Q-Marine streamers under the 2.5D earth model, thefilter has three coefficients only. Thus, to estimate the verticalparticle velocity at a particular receiver location, only the pressurefield at this and the two neighbouring receivers are required. Providedthe depth of the (near-horizontal) streamers below the sea surface areknown at all points along the streamers, equation (A1) gives amethodology of estimating the vertical particle velocity field inrough-sea environments. Let P and V_(z) denote the spatial Fouriertransforms of p and v_(z), respectively. The estimate given in equation(A1) then can be combined with the raw pressure measurements to deghostthe data (see e.g. L. Amundsen in Geophysics Vol. 56, pp 1027-1039(1993) or UK patent application No. 9906456.0) $\begin{matrix}{{{\overset{\sim}{P}\left( {\omega,k_{x},k_{y},z} \right)} = {\frac{1}{2}\left\lbrack {{P\left( {\omega,k_{x},k_{y},z} \right)} - {\frac{\rho\quad\omega}{k_{z}}{V_{z}\left( {\omega,k_{x},k_{y},z} \right)}}} \right\rbrack}},} & ({A3})\end{matrix}$where {tilde over (P)}(ω,k_(x),k_(y),z) denotes the deghosted (up-going)pressure, k_(z)=√{square root over (k²−κ² )} is the vertical wavenumber,κ²=k_(x) ²+k_(y) ² and k_(x) and k_(y) are the horizontal wavenumbers.Equation (3) is valid for all wavenumbers. Alternatively, approximationsto equation (A3) as those suggested by A. Osen et al. in “Towardsoptimal spatial filters for demultiple and wavefield splitting of oceanbottom seismic data”, submitted to Geophysics (2002).

In the present application the accuracy of the vertical particlevelocity field approximation given in equation (A1) is extended byintroducing optimised filters. L. Amundsen et al. derived, in GeophysicsVol. 60, pp 212-222 (1995) an exact frequency-wavenumber domainrelationship between the vertical component of the particle velocity andthe pressure field. For single streamer data the relationship is:$\begin{matrix}{{V_{z}\left( {\omega,k_{x},k_{y},z} \right)} = {{- \frac{k_{z}}{\rho\quad\omega}}{\left( \frac{{\exp\left( {{- {ik}_{z}}z} \right)} + {\exp\left( {{ik}_{z}z} \right)}}{{\exp\left( {{- {ik}_{z}}z} \right)} - {\exp\left( {{ik}_{z}z} \right)}} \right) \cdot {P\left( {\omega,k_{x},k_{y},z} \right)}}}} & ({A4})\end{matrix}$

Equation (A4) is valid for all streamer depths z in the water columnwhen the receiver cable is horizontal and the air/water surface is flatwith vanishing pressure. Essentially, the equation shows that to findthe vertical component of the particle velocity from the pressure field,first the pressure has to be receiver-deghosted, then the receiver ghostoperator for the vertical component of the particle velocity is applied.(Taking into account that the reflection coefficient of the free surfaceis −1, the ghost operator for the pressure is 1−exp(2ik_(z)z) since thepressure is the sum of upgoing and downgoing waves, whereas the ghostoperator for the vertical particle velocity is 1+exp(2ik_(z)z) since thevertical particle velocity is proportional to the difference betweenupgoing and downgoing waves.) Note that when the source depth is lessthan the receiver depth the incident pressure wavefield (the sum of thedirect pressure wavefield and its ghost) does not contain a receiverghost, but only a source ghost. Since equation (A4) relies on receiverdeghosting, it cannot correctly treat the incident wavefield. This is ofminor importance in receiver deghosting of streamer data. When thesource depth is greater than the receiver depth, however, the incidentpressure wavefield contains a receiver ghost as do the reflected part ofthe pressure field. In this case, equation (A4) is valid for the fallpressure wavefield.

Equation (A4) can be rewritten as $\begin{matrix}{{{V_{z}\left( \quad{\omega,\quad k_{x},\quad k_{y},\quad z} \right)} = \quad{{- \frac{i}{\rho\quad\omega\quad z}}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}\quad{F_{m}\left( \quad{k,\quad z} \right)}\quad{\kappa^{2m} \cdot \quad{P\left( \quad{\omega,\quad k_{x},\quad k_{y},\quad z}\quad \right)}}}}}}\quad,} & ({A5})\end{matrix}$with horizontal wavenumber (κ) independent functions $\begin{matrix}{F_{m} = {{F_{m}\left( {k,z} \right)} = {k^{{- 2}m}{\sum\limits_{n = m}^{\infty}{A_{nm}({kz})}^{2n}}}}} & ({A6}) \\{where} & \quad \\{A_{nm} = {\frac{\left( {- 1} \right)^{n}2^{2n}B_{2n}}{\left( {2n} \right)!}\begin{pmatrix}n \\m\end{pmatrix}}} & ({A7})\end{matrix}$are scalars, B_(n) is the n:th Bernoulli number, and n! is the factorialof n. All functions F_(m) can be given in compact form. For instance,$\begin{matrix}{F_{0} = {{{kz}\quad{\cot({kz})}} = {{- {\mathbb{i}}}\quad{{kz}\left( \frac{{\exp\left( {- {ikz}} \right)} + {\exp({ikz})}}{{\exp\left( {- {ikz}} \right)} - {\exp({ikz})}} \right)}}}} & ({A8}) \\{{F_{1} = {\frac{z}{2k}\frac{\partial}{\partial({kz})}F_{0}}},} & ({A9}) \\{F_{2} = {\frac{z}{4k^{3}}\left( {{\frac{\partial}{\partial({kz})}{kz}} - 2} \right)\quad\frac{\partial}{\partial({kz})}{F_{0}.}}} & ({A10})\end{matrix}$

The function F₀ ensures that vertically propagating waves(k_(x)=k_(z)=0) are correctly deghosted.

Taking an inverse Fourier transform of equation (A5) over horizontalwavenumbers gives $\begin{matrix}{{{v_{z}\left( {\omega,x,y,z} \right)} = {{- \frac{\mathbb{i}}{{\rho\omega}\quad z}}{\sum\limits_{m = 0}^{\infty}{{F_{m}\left( {k,z} \right)}\quad{D_{2}^{m}\left( {x,y} \right)}^{*}{p\left( {\omega,x,y,z} \right)}}}}},} & ({A11})\end{matrix}$where D₂ is the inverse Fourier transform of κ², and * is the 2D spatialconvolution operator. For infinite bandwidth, $\begin{matrix}{{{D_{2}\left( {x,y} \right)} = {{D_{2}(x)} + {D_{2}(y)}}},} & ({A12}) \\{where} & \quad \\\begin{matrix}{{{D_{2}(x)} = \frac{\partial^{2}}{\partial x^{2}}};} & {{{D_{2}(y)} = \frac{\partial^{2}}{\partial y^{2}}},}\end{matrix} & ({A13})\end{matrix}$are second-order spatial differentiation operators. Differentiationhowever involves large amplification of high wavenumbers. Hence, inpractical situations a bandlimited version of differentiation ispreferably used. Let d₂ denote a bandlimited version of D₂. Two possiblebandlimited differentiation operators are $\begin{matrix}{{{d_{2}(x)} = {\frac{1}{\Delta\quad x}\left\lbrack {{\delta\left( {x + {\Delta\quad x}} \right)} - {2\quad{\delta(x)}} + {\delta\left( {x - {\Delta\quad x}} \right)}} \right\rbrack}},} & ({A14})\end{matrix}$where Δx is the spatial sampling interval, or $\begin{matrix}{{{d_{2}\left( {x,\omega} \right)} = {\frac{1}{\pi}{\int_{0}^{K}{{\mathbb{d}k_{x}}\quad{\cos\left( {k_{x}x} \right)}\quad{W\left( {\omega,k_{x}} \right)}\left( {- k_{x}^{2}} \right)}}}},} & ({A15})\end{matrix}$where W is a weighting function that properly bandlimits thedifferentiation operator. Note that W can be changed per frequencycomponent. In equation (A15) K is the maximum wavenumber. Two possiblechoices are K=ω/c, where c is the velocity of water, or K=π/Δx, theNyquist wavenumber.

Note that equation (A11) is valid for all receiver depths z in the watercolumn. Also, observe that the scheme is iterative:D ₂ ^(m) *p=D ₂ ^(m)*(D ₂ ^(m−1) *p);D ₂ ⁰=1  (A16)

This recursive property allows equation (A11) to be implemented withnumerical efficiency. For dense spatial sampling the differentiatorsD₂(x) and D₂(y) can be approximated by filters of length three [seeequation (A14)]. A properly designed, bandlimited, second-orderdifferentiator d₂(x) is imperative for obtaining a stable iterativescheme.

Special Case I: Finite m

A special case of equation (A11) arises when m is finite, say M. Then,for infinite bandwidth, $\begin{matrix}{{{v_{z}\left( {\omega,x,y,z} \right)} \approx {{- \frac{\mathbb{i}}{{\rho\omega}\quad z}}{\sum\limits_{m = 0}^{M}{{F_{m}\left( {k,z} \right)}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)^{m}\quad{p\left( {\omega,x,y,z} \right)}}}}},} & ({A17})\end{matrix}$ where F_(m) is given in equation (A6) as an infiniteseries. However, the function F_(m) can also be given analytically interms of trigonometric functions. Examples are F₀, F₁, and F₂ given inequations (A8), A(9), and (A10), respectively. Observe that theassumption of flat sea surface used to derive equation (4) can now berelaxed. In equation (A17) the x and y variables are decoupled from z,implying that it is valid also for a rough sea surface.

Equation (A17) can be numerically implemented in different ways. Thepreferred implementation depends on the spatial sampling interval forthe pressure field. For densely sampled data one would typicallyrepresent the spatial derivatives by classical central finite-differenceapproximations. If M is small, the filter is then local and compact. Asan example, choosing M=2 gives for infinite bandwidth, $\begin{matrix}{{v_{z}\left( {\omega,x,y,z} \right)} \approx {{- {\frac{\mathbb{i}}{{\rho\omega}\quad z}\left\lbrack {F_{0} + {F_{1}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)} + {F_{2}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)}^{2}} \right\rbrack}}\quad{{p\left( {\omega,x,y,z} \right)}.}}} & ({A18})\end{matrix}$

Approximating the second-order spatial derivatives by a three-pointfinite difference operator the iterative filter is of length three.Since the functions F₀, F₁and F₂ are valid for all receiver depths,equation (A18) is more general than equation (A1), and is preferred forrough-sea deghosting.

For a coarser receiver interval the horizontal derivatives are mostconveniently implemented as numerically optimised difference operatorsor by use of Fourier techniques.

Special Case II: Finite m, n

Limiting the sums over m in equation (11) and it in equation (A6) frominfinity to M and N(m), respectively, gives for infinite bandwidth theapproximation with $\begin{matrix}{{{v_{z}\left( {\omega,x,y,z} \right)} \approx {{- \frac{\mathbb{i}}{\rho\quad\omega\quad z}}{\sum\limits_{m = 0}^{M}{{F_{m}\left( {k,z} \right)}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)^{m}{p\left( {\omega,x,y,z} \right)}}}}},} & ({A19}) \\{with} & \quad \\{F_{m} \approx {k^{{- 2}m}{\sum\limits_{n = m}^{N{(m)}}{A_{nm}({kz})}^{2n}}}} & ({A20}) \\{{{{Choosing}\quad M} = 1},{{N(0)} = 2},{{{and}\quad{N(1)}} = 1},{gives}} & \quad \\{{F_{0} \approx {1 - \frac{({kz})^{2}}{3}}},} & ({A21}) \\{F_{1} \approx {- \frac{z^{2}}{3}}} & ({A22}) \\{and} & \quad \\{{{v_{z}\left( {\omega,x,y,z} \right)} \approx {{- {\frac{\mathbb{i}}{\rho\quad\omega\quad z}\left\lbrack {A_{00} + {A_{10}({kz})}^{2} + {A_{11}{z^{2}\left( {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}} \right)}}} \right\rbrack}}\quad{p\left( {\omega,x,y,z} \right)}}},} & ({A23}) \\{{with}\quad{scalars}} & \quad \\\begin{matrix}{{A_{00} = 1};} & {{A_{10} = {{- 1}/3}};} & {A_{11} = {{- 1}/3.}}\end{matrix} & (24)\end{matrix}$

Equation (A23) has the same structure as equation (A1) proposed byRobertsson and Kragh (above). The scalars A₁₀ and A₁₁, however, differfrom the scalars A₁₀ ^({acute over ( )}) and A₁₁ ^({acute over ( )}) inequation (A1). In many respects the limitations of equation (A23) aresimilar to those of the formula of Robertsson and Kragh (above).However, equations (A11), (A17), and (A19) are general and onesignificant difference is that they are valid for all streamer depths,whereas the theory derived by Robertsson and Kragh (above) requires thatthe streamer is towed no deeper than approximately λ/3.5, where λ is theminimum wavelength.

If we choose M=1, N(0)=∞ and N(1)=2, this results in F₀=kz cot(kz) and$F_{1} = {{- \frac{z^{2}}{3}}{\left( {1 + {\frac{2}{15}k^{2}z^{2}}} \right).}}$Data Processing of Single Streamer Data

Equations (A11), (A17) and (A19) above include derivatives in they-direction or cross-line direction (the streamer(s) is/are assumed toextend along the x-direction). In a seismic surveying arrangement thatincludes only a single streamer it is not clear how the cross-linederivatives in the y-direction that occur in the above equations may bedetermined from the seismic data acquired by the receivers on thestreamer. It will be noted, however, that only even powers(second-order, fourth-order, sixth-order, etc.) spatial derivativesoccur. Finding an expression for the second-order spatial derivative istherefore sufficient.

Consider the single streamer geometry shown in FIG. 4, which is aschematic plan view of the seismic surveying arrangement of FIG. 1. Theseismic source 1 is located at coordinates (O,y₀) in the x-y plane. Thestreamer extends along the x-axis, and has a y-coordinate y₀. One of thereceivers S_(i) is located at coordinates (x₀,y₀), and it assumed thatthe source 1 and the receiver S_(i) are located in the same z-plane inan acoustic layer. It will be assume that the model only varies withdepth, and does not vary with x or y. By using a finite-differenceapproximation to a second-order derivative, it is possible to write thefollowing equation for the pressure p: $\begin{matrix}{{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {{\frac{1}{\Delta\quad y^{2}}\left( {{p\left( {x_{0},{y_{0} + {\Delta\quad y}}} \right)} - {2\quad{p\left( {x_{0},y_{0}} \right)}} + {p\left( {x_{0},{y_{0} - {\Delta\quad y}}} \right)}} \right)} + {O\left( {\Delta\quad y^{2}} \right)}}},} & ({A25})\end{matrix}$where O(Δy²) expresses the leading error term in the finite-differenceapproximation. However, because of radial symmetry the recorded pressureis the same on either side of the receiver so that equation (A25)becomes $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {{\frac{2}{\Delta\quad y^{2}}\left( {{p\left( {x_{0},{y_{0} + {\Delta\quad y}}} \right)} - {p\left( {x_{0},y_{0}} \right)}} \right)} + {{O\left( {\Delta\quad y^{2}} \right)}.}}} & ({A26})\end{matrix}$

Radial symmetry also yieldsp(x ₀ ,y ₀ +Δy)=p(x ₀√{square root over (1+(Δy/x ₀)²)},y ₀)  (A27)so that equation (A26) becomes: $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {{\frac{2}{\Delta\quad y^{2}}\left( {{p\left( {{x_{0}\sqrt{1 + \left( {\Delta\quad{y/x_{0}}} \right)^{2}}},y_{0}} \right)} - {p\left( {x_{0},y_{0}} \right)}} \right)} + {{O\left( {\Delta\quad y^{2}} \right)}.}}} & ({A28})\end{matrix}$

By using a first-order one-sided finite-difference approximation in thex -direction it is possible to write equation (A28) as: $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {{\frac{2{x_{0}\left( {\sqrt{1 + \left( {\Delta\quad{y/x_{0}}} \right)^{2}} - 1} \right)}}{\Delta\quad y^{2}}\left( {{\frac{\partial}{\partial x}{p\left( {x_{0},y_{0}} \right)}} + {O\left( {\sqrt{1 + \left( {\Delta\quad{y/x_{0}}} \right)^{2}} - 1} \right)}} \right)} + {O\left( {\Delta\quad y^{2}} \right)}}} & ({A29})\end{matrix}$

Taking a Taylor series expansion of the square root: $\begin{matrix}{{\sqrt{1 + \left( {\Delta\quad{y/x_{0}}} \right)^{2}} = {1 + {\frac{1}{2}\left( {\Delta\quad{y/x_{0}}} \right)^{2}} + {O\left( {\Delta\quad y^{4}} \right)}}},} & ({A30})\end{matrix}$equation (A29) may be re-written as, $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {{\frac{2\quad{x_{0}\left( {{\frac{1}{2}\left( {\Delta\quad{y/x_{0}}} \right)^{2}} + {O\left( {\Delta\quad y^{4}} \right)}} \right)}}{\Delta\quad y^{2}}\left( {{\frac{\partial}{\partial x}{p\left( {x_{0},y_{0}} \right)}} + {O\left( {\Delta\quad y^{2}} \right)}} \right)} + {{O\left( {\Delta\quad y^{2}} \right)}.}}} & ({A31})\end{matrix}$

Further simplification of equation (31) yields $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {{\frac{1}{x_{0}}\frac{\partial}{\partial x}{p\left( {x_{0},y_{0}} \right)}} + {{O\left( {\Delta\quad y^{2}} \right)}.}}} & ({A32})\end{matrix}$

Finally, we let Δy→0 so that: $\begin{matrix}{{\frac{\partial^{2}}{\partial y^{2}}{p\left( {x_{0},y_{0}} \right)}} = {\frac{1}{x_{0}}\frac{\partial}{\partial x}{{p\left( {x_{0},y_{0}} \right)}.}}} & ({A33})\end{matrix}$

For pressure data over a horizontally layered medium there is no problemwith the singularity in equation (33) since$\frac{\partial}{\partial x}{p\left( {x_{0},y_{0}} \right)}$goes to zero faster than the singularity 1/x₀ goes to infinity.

Equation (A33) is consistent with what is valid for cylindrical symmetry(data over a horizontally layered medium). An alternative approach wouldtherefore be to multiply the data with √{square root over (t)}, where tis time, and use the 2D version of the theory derived in thisapplication as well as in UK patent application No. 0015810.5. However,equation (A33) can potentially give more accurate results.

It should be noted the above approach may not be reliable for a mediumthat is not simply plane layered since there can be a problem because pis not symmetric with respect to x.

Complex Frequency

In the implementation of the filters we use complex frequencies to avoidthe singularities in the functions F_(m) as given e.g. in equations (8),(9), and (10) for M=0,1,2, respectively. The use of a complex-frequencyis described by, for example, S. Mallick and L. N. Frazer in GeophysicsVol. 52, pp 1355-1364 (1987) and by L. Amundsen and B. Ursin inGeophysics Vol. 56, pp 1027-1039 (1991).

On V_(z)-Estimation using Numerically Optimised Spatial Filters

1. Introduction

In this section, we introduce numerically optimised spatial filters forthe estimation of the vertical particle velocity field. Amundsen et al.(1995, above) derived an exact frequency-wavenumber domain relationshipbetween the vertical component of the particle velocity and the pressurefield. For single streamer data the relationship in the 2D case is$\begin{matrix}{{v_{z}\left( {\omega,k_{x},z} \right)} = {{- {\frac{k_{z}}{\rho\quad\omega}\left\lbrack \frac{{\exp\left( {{- {ik}_{z}}z} \right)} + {\exp\left( {{ik}_{z}z} \right)}}{{\exp\left( {{- {ik}_{z}}z} \right)} - {\exp\left( {{ik}_{z}z} \right)}} \right\rbrack}} \cdot {p\left( {\omega,k_{x},z} \right)}}} & ({A34})\end{matrix}$

Equation (A34) may be re-written as:v _(z)(ω, k _(x) , z)=F(ω, k _(x) ,z)·p(ω, k _(x) , z)  (A35)where the v_(z)-estimation filter is given by: $\begin{matrix}{{F\left( {\omega,k_{x},z} \right)} = {- {\frac{k_{z}}{\rho\quad\omega}\left\lbrack \frac{{\exp\left( {{- {ik}_{z}}z} \right)} + {\exp\left( {{ik}_{z}z} \right)}}{{\exp\left( {{- {ik}_{z}}z} \right)} - {\exp\left( {{ik}_{z}z} \right)}} \right\rbrack}}} & ({A36})\end{matrix}$

In the frequency-space domain, equation (A34) is written symbolically asV _(z) ={tilde over (F)}*p  (A37)where * denotes spatial convolution.2. Optimisation Problem FormulationThe discrete horizontal wavenumber response of a spatial digital filtercan be written as: $\begin{matrix}{{\overset{\sim}{F}}_{k} = {{\overset{\sim}{F}\left( {k\quad\Delta\quad k_{x}} \right)} = {\frac{\sum\limits_{m = {{- M}/2}}^{M/2}{b_{m}\quad{\exp\left( {{- {\mathbb{i}}}\quad k\quad\Delta\quad k_{x}m\quad\Delta\quad x} \right)}}}{\sum\limits_{n = 0}^{N}{a_{n}\quad{\exp\left( {{- {\mathbb{i}}}\quad k\quad\Delta\quad k_{x}n\quad\Delta\quad x} \right)}}} = \frac{\Theta_{km}b_{m}}{\Phi_{kn}a_{n}}}}} & ({A38})\end{matrix}$where a_(n),b_(m) are the backward and forward filter coefficients withN,M the respective filter orders. Without loss of generality, M and Nare assumed to be even numbers and a₀≡1. Furthermore, Δx is the spatialsampling interval and Δk_(x)=π/[Δx(K−1)] is the distance between thediscrete horizontal wavenumbers. If F_(k)=F(kΔk_(x)) denotes the desireddiscrete horizontal wavenumber response, then the design of the spatialfilters can be stated as a general weighted non-linear least-squaresoptimisation problem: $\begin{matrix}{{\overset{\sim}{F}}_{k} = {\min\limits_{b_{m},a_{n}}{\sum\limits_{k = {- {({K - 1})}}}^{K - 1}{W_{k}\quad{{\frac{\Theta_{km}b_{m}}{\Phi_{kn}a_{n}} - F_{k}}}^{2}}}}} & ({A39})\end{matrix}$where W_(k) is a set of positive weights. In the following, theoptimisation problem is formulated for propagating waves only, for whichthe ideal decomposition filter is purely real (i.e., zero-phase). Let κdenote the index corresponding to a discrete horizontal wavenumber closeto the critical wavenumber.3. Zero-Phase FIR Filters

For a zero-phase non-recursive or finite impulse response (FIR) filter,the forward filter coefficients are symmetric and the backward filterorder is zero. This reduces the non-linear least-squares optimisationproblem given in equation (A39) to the following linear-least squaresoptimisation problem: $\begin{matrix}{\min\limits_{g_{m}}{\sum\limits_{k = 0}^{\kappa - 1}{W_{k}{{{\Gamma_{km}g_{m}} - F_{k}}}^{2}}}} & ({A40})\end{matrix}$where Γ_(km)=cos(kΔk_(x)mΔx). The g_(m)'s are related to the b_(m)'s byb_(m)=0.5 g_(m) for g_(m)=1,2, . . . , M/2 and by b₀=g₀.4. Quadratic Programming

The linear least-squares optimisation problem shown in equation (A40) isunconstrained, and unique and optimal solutions are given by thepseudo-inverses. Nevertheless, it can be important that the magnituderesponses of the optimal FIR filter have zero-error at particularhorizontal wavenumbers, or that they have magnitude responses smallerthan a given function for a certain horizontal wavenumber range. Thecombination of a linear least-squares optimisation problem with linearequality and/or inequality constraints is called Quadratic Programming(QP) and can be solved efficiently with many different algorithms.

For example, a single equality constraint can be included at zerohorizontal wavenumber to ensure that propagating waves are estimatedcorrectly by the optimal FIR filter. In addition, linear inequalityconstraints can be included in the evanescent horizontal wavenumberrange to force the spatial filter's magnitude response of the idealdecomposition filter. This prevents boosting of evanescent waves by thenumerically optimised spatial filters.

5. Twin Streamers

In this section we describe how to filter data recorded in twin-streamerconfigurations to obtain an estimate of particle velocity at theuppermost streamer. The advantage of the proposed technique is that acompact filter is achieved which is suitable when the streamerseparation varies along the streamers. Once the particle velocity hasbeen estimated, it can be combined with the pressure recording asdescribed by for instance PCT application PCT/GB00/01074 (see equation(A43) below) to separate the recorded data into its up- and down-goingcomponents.

Theory

L. Amundsen gives, at “Wavenumber-based filtering of marine point-sourcedata”, Geophysics, Vol. 58, pp 1335-1348 (1993), an expression for theup-going pressure {tilde over (p)}(ω, k_(x), k_(y), z₁) at theshallowest streamer in a twin-streamer configuration with two horizontalstreamers with a constant vertical separation: $\begin{matrix}{{{\overset{\sim}{p}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = \frac{{{\exp\left( {{\mathbb{i}}\quad k_{z}\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}} - {{\exp\left( {2\quad{\mathbb{i}}\quad k_{z}\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}}{1 - {\exp\left( {2{\mathbb{i}}\quad k_{z}\Delta\quad z} \right)}}},} & ({A41})\end{matrix}$where p(ω, k_(x), k_(y), z₂ ) and p(ω, k_(x), k_(y), z₁) denote therecorded pressure at the deepest and shallowest streamers, respectively,Δz is the vertical streamer separation, ω is the angular frequency,k_(z)=√{square root over (k²−κ²)} is the vertical wavenumber, κ²=k_(x)²+k_(y) ², and k_(x) and k_(y) are the horizontal wavenumbers. Equation(A41) can be expanded as: $\begin{matrix}{{\overset{\sim}{p}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {\frac{1}{2}{\left( {{p\left( {\omega,k_{x},k_{y},z_{1}} \right)} + {\frac{\mathbb{i}}{\sin\left( {k_{z}\Delta\quad z} \right)}{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}} - {{\mathbb{i}}\quad\cot\quad\left( {k_{z}\Delta\quad z} \right)\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}} \right).}}} & ({A42})\end{matrix}$

The up-going pressure at the shallowest streamer can also be expressedin terms of pressure and vertical particle velocity recordings v_(z)(ω,k_(x), k_(y), z₁) at the shallowest streamer as: $\begin{matrix}{{{\overset{\sim}{p}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {\frac{1}{2}\left\lbrack {{p\left( {\omega,k_{x},k_{y},z_{1}} \right)} - {\frac{\rho\quad\omega}{k_{z}}{v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)}}} \right\rbrack}},} & ({A43})\end{matrix}$(e.g., Amundsen 1993 (above) or PCT/GB00/01074) where ρ is the density.By identification of terms between equations (A42) and (A43) we obtainan expression for vertical particle velocity along the shalloweststreamer in terms of recorded pressure along the two streamers in a twinstreamer configuration: $\begin{matrix}{{v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{{\mathbb{i}}\quad k_{z}}{\rho\omega}\quad{\cot\left( {k_{z}\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}} - {\frac{{\mathbb{i}}\quad k_{z}}{\rho\quad\omega\quad{\sin\left( {k_{z}\Delta\quad z} \right)}}\quad{{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}.}}}} & ({A44})\end{matrix}$

Equation (A44) closely resembles a corresponding expression for singlestreamer de-ghosting obtained by L. Amundsen et al.: $\begin{matrix}{{v_{z}\left( {\omega,k_{x},k_{y},z} \right)} = {{- \frac{{\mathbb{i}}\quad k_{z}}{\rho\quad\omega}}{\cot\left( {k_{z}\Delta\quad z} \right)}{{p\left( {\omega,k_{x},k_{y},z} \right)}.}}} & ({A45})\end{matrix}$

Therefore, all that remains is to establish compact filters for thesecond term in equation (A44), which we denote as: $\begin{matrix}{G^{(2)} = {- {\frac{{\mathbb{i}}\quad k_{z}}{\rho\quad\omega\quad\sin\quad\left( {k_{z}\Delta\quad z} \right)}.}}} & ({A46})\end{matrix}$

According to Raade and Westergren in Beta—Mathematics handbook:Studentlitteratur, Lund, Sweden (1988), equation (A46) can be expandedas: $\begin{matrix}{{G^{(2)} = {{- \frac{\mathbb{i}}{{\rho\omega\Delta}\quad z}}{\sum\limits_{n = 0}^{\infty}{\left( {- 1} \right)^{n - 1}\frac{2^{2n} - 2}{\left( {2n} \right)!}{B_{2n}\left( {k_{z}\Delta\quad z} \right)}^{2n}}}}},} & ({A47})\end{matrix}$where B_(n) is the nth Bernoulli number and −π<k_(z)Δz<π. Making use ofthe Binomial expansion to expand factors (k_(z)Δz)² as outlined inAmundsen et al. (2003), equation $\begin{matrix}{{G^{(2)} = {{- \frac{\mathbb{i}}{{\rho\omega\Delta}\quad z}}{\sum\limits_{n = 0}^{\infty}{\left( {- 1} \right)^{n - 1}\frac{2^{2n} - 2}{\left( {2n} \right)!}{B_{2n}\left( {k\quad\Delta\quad z} \right)}^{2n}\quad{\sum\limits_{m = 0}^{n}{\left( {- 1} \right)^{m}\begin{pmatrix}n \\m\end{pmatrix}\quad k^{{- 2}m}\kappa^{2m}}}}}}},} & ({A48})\end{matrix}$where k=ω/c and c is the speed of sound in water. Interchanging sumsusing (using $\begin{matrix}{\left. {{\sum\limits_{n = 0}^{\infty}{a_{n}\quad{\sum\limits_{m = 0}^{n}{b_{nm}c_{m}}}}} = {\sum\limits_{m = 0}^{\infty}{c_{m}\quad{\sum\limits_{n = m}^{\infty}{a_{n}b_{nm}}}}}} \right)\quad{yields}\text{:}} & \quad \\{{G^{(2)} = {{- \frac{i}{\rho\quad\omega\quad\Delta\quad z}}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}k^{{- 2}m}\kappa^{2m}\quad{\sum\limits_{n = m}^{\infty}{A_{nm}\left( {k\quad\Delta\quad z} \right)}^{2n}}}}}},} & ({A49}) \\{where} & \quad \\{A_{nm} = {\left( {- 1} \right)^{n - 1}\frac{2^{2n} - 2}{\left( {2n} \right)!}{{B_{2n}\begin{pmatrix}n \\m\end{pmatrix}}.}}} & ({A50})\end{matrix}$

By introducing$F_{m} = {k^{{- 2}m}{\sum\limits_{n = m}^{\infty}{A_{nm}\left( {k\quad\Delta\quad z} \right)}^{2n}}}$we can rewrite equation (A44) as $\begin{matrix}{{v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{i}{\rho\quad\omega\quad\Delta\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(1)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}}} - {\frac{i}{{\rho\omega\Delta}\quad z}{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(2)}\kappa^{2m}{{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}.}}}}}} & ({A51})\end{matrix}$

The lowest order coefficients in the first expansion are given inAmundsen et al. (2003): $\begin{matrix}{{F_{0}^{(1)} = {k\quad\Delta\quad z\quad{\cot\left( {k\quad\Delta\quad z} \right)}}},} & ({A52}) \\{and} & \quad \\{F_{1}^{(1)} = {\frac{\Delta\quad z}{2k}\frac{\mathbb{d}}{\mathbb{d}\left( {k\quad\Delta\quad z} \right)}k\quad\Delta\quad z\quad{{\cot\left( {k\quad\Delta\quad z} \right)}.}}} & ({A53})\end{matrix}$

The two lowest order coefficients in the second expansion can beidentified to be: $\begin{matrix}{{F_{0}^{(2)} = \frac{k\quad\Delta\quad z}{\sin\left( {k\quad\Delta\quad z} \right)}},} & ({A54}) \\{and} & \quad \\{F_{1}^{(2)} = {\frac{\Delta\quad z}{k}\quad{\frac{{\sin\left( {k\quad\Delta\quad z} \right)} - {k\quad\Delta\quad z\quad{\cos\left( {k\quad\Delta\quad z} \right)}}}{\sin^{2}\left( {k\quad\Delta\quad z} \right)}.}}} & ({A55})\end{matrix}$

By only keeping the lowest order terms [equations (A52) and (A54)] inthe series expansions in equation (A51) results in an expression that isvalid for all streamer separations and that is exact for verticallyincident waves only. This is equivalent to the known as the “shift &sum” technique used in the early days of twin-streamer deghosting.

The results can be significantly improved by keeping higher order termsin the series expansions. The next higher order approximation willinclude 4 terms as given by equations (A52)-(A55). This can beimplemented with a 3-point running spatial filter.

1. A method of determining a vertical component of particle velocityfrom seismic data acquired using at least one receiver disposed within awater column, the method comprising: determining the vertical componentof particle velocity from the pressure at a first location and a secondlocation, the first location being the location of a first receiver andbeing vertically beneath the second location, using an operator suitablefor seismic data acquired at a vertical separation between the first andsecond locations of up to at least 0.4 times the minimum wavelength ofseismic energy emitted by a seismic source used to acquire the seismicdata.
 2. A method as claimed in claim 1 wherein the second location isat or above the surface of the water column.
 3. A method as claimed inclaim 1 wherein the second location is below the surface of the watercolumn and is the location of a second receiver.
 4. A method as claimedin claim 3 and comprising determining the vertical component of particlevelocity at the second location from pressure acquired at the secondreceiver and from pressure acquired at the first receiver.
 5. A methodas claimed in claim 4 and comprising determining the vertical componentof particle velocity at the second location according to the followingequation or from an equation derived therefrom:${v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{{ik}_{z}}{\rho\quad\omega}{\cot\left( {k_{z}\quad\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}} - {\frac{{ik}_{z}}{\rho\quad\omega\quad\sin\quad\left( {k_{z}\Delta\quad z} \right)}{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}}}$where z₁ denotes the depth below the surface of the water column of thesecond location, z₂ denotes the depth of the first location, Δz=z₂−z₁, pdenotes the pressure, v_(z) denotes vertical component of particlevelocity, ρ denotes the density of the water column, ω denotes theangular frequency, k_(x) and k_(y) denote the horizontal wavenumber andk_(z) denotes the vertical wavenumber.
 6. A method as claimed in claim 4or 5 and comprising determining the vertical component of particlevelocity at the second location from the following equation:${v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{i}{{\rho\omega\Delta}\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}\quad F_{m}^{(1)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}}} - {\frac{i}{{\rho\omega\Delta}\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(2)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}}}}}$$\begin{matrix}{where} & {{F_{m} = {k^{{- 2}m}\quad{\sum\limits_{n = m}^{\infty}{A_{nm}\left( {k\quad\Delta\quad z} \right)}^{2n}}}},} & {{A_{nm} = {\left( {- 1} \right)^{n - 1}\frac{2^{2n} - 2}{\left( {2n} \right)!}{B_{2n}\begin{pmatrix}n \\m\end{pmatrix}}}},}\end{matrix}$ B_(n) is the n^(th) Bernoulli number, and κ²=k_(x) ²+k_(y)².
 7. A method as claimed in claim 6 wherein each summation over m iscarried out over m=0 to m=1.
 8. A method as claimed in any of claim 1 to7 wherein the operator is a spatially compact operator.
 9. A method asclaimed in any of claims 1 to 8 and comprising the further step ofprocessing the seismic data using the vertical component of the particlevelocity.
 10. A method as claimed in claim 9 and comprising processingthe seismic data so as to determine at least one of an up-goingcomponent and a down-going component of the seismic data.
 11. A methodas claimed in claim 9 and comprising processing the seismic data usingthe vertical component of the particle velocity thereby to attenuateeffects in the processed data of seismic energy reflected and/orscattered at the free-surface of the water column.
 12. A method asclaimed in claim 9 and comprising processing the seismic data using thevertical component of the particle velocity thereby to obtaininformation about the signature of the source of seismic energy.
 13. Amethod as claimed in any of claims 1 to 8 wherein the seismic datacomprise pressure and particle velocity data, and wherein the methodfurther comprises comparing the vertical component of the particlevelocity determined from the acquired pressure with the measured valuesof the particle velocity.
 14. A method as claimed in any of claims 1 to8 wherein the seismic data comprises pressure and particle velocitydata, and wherein the method further comprises the step of determiningthe depth of the seismic receiver within the water column from themeasured vertical component of the particle velocity and the verticalcomponent of the particle velocity determined from the pressure.
 15. Amethod of determining a vertical component of particle velocity fromseismic data acquired at first and second vertically-separated streamersdisposed within a water column, the first streamer comprising M (M>1)pressure sensors disposed along the length of the first streamer and thesecond streamer comprising N (N>1) pressure sensors disposed along thelength of the second streamer, the method comprising applying to thepressure acquired at sensors on at least one streamer an operator havinga length (in the spatial domain) that is greater than one but that isless than the number of pressure sensors on the respective streamer. 16.A method of acquiring marine seismic data comprising the steps of:actuating an array of one or more seismic sources to emit seismicenergy; acquiring seismic data at one or more receivers disposed withina water column, the seismic data including at least pressure data; andprocessing the acquired pressure data according to a method as definedin any of claims 1 to
 15. 17. An apparatus for determining a verticalcomponent of particle velocity from seismic data acquired using at leastone receiver disposed within a water column, the apparatus having meansfor determining the vertical component of particle velocity from thepressure at a first location and a second location, the first locationbeing the location of a first receiver and being vertically beneath thesecond location, using an operator suitable for seismic data acquired ata vertical separation between the first and second locations of up to atleast 0.4 times the minimum wavelength of seismic energy emitted by aseismic source used to acquire the seismic data.
 18. An apparatus asclaimed in claim 17 and adapted to determine the vertical component ofparticle velocity from the pressure acquired at the first and secondlocations, the second location being at or above the surface of thewater column.
 19. An apparatus as claimed in claim 17 and adapted todetermine the vertical component of particle velocity from the pressureacquired at the first and second locations, the second location beingwithin the water column and being the location of a second receiver. 20.An apparatus as claimed in claim 19 and adapted to determine thevertical component of particle velocity at the second location frompressure acquired at the second receiver and from pressure acquired atthe first receiver.
 21. An apparatus as claimed in claim 20 and adaptedto determine the vertical component of particle velocity at the secondlocation according to the following equation or from an equation derivedtherefrom:${v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{{ik}_{z}}{\rho\quad\omega}{\cot\left( {k_{z}\quad\Delta\quad z} \right)}\quad{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}} - {\frac{{ik}_{z}}{\rho\quad\omega\quad\sin\quad\left( {k_{z}\Delta\quad z} \right)}{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}}}$where z₁ denotes the depth below the surface of the water column of thesecond location, z₂ denotes the depth of the first location, Δz=z₂−Z₇, pdenotes the pressure, v_(z) denotes vertical component of particlevelocity, ρ denotes the density of the water column, ω denotes theangular frequency, k_(x) and k_(y) denote the horizontal wavenumber andk_(z) denotes the vertical wavenumber.
 22. An apparatus as claimed inclaim 20 or 21 and adapted to determine the vertical component ofparticle velocity at the second location from the following equation:${v_{z}\left( {\omega,k_{x},k_{y},z_{1}} \right)} = {{\frac{i}{{\rho\omega\Delta}\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}\quad F_{m}^{(1)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{1}} \right)}}}} - {\frac{i}{{\rho\omega\Delta}\quad z}\quad{\sum\limits_{m = 0}^{\infty}{\left( {- 1} \right)^{m}F_{m}^{(2)}\kappa^{2m}{p\left( {\omega,k_{x},k_{y},z_{2}} \right)}}}}}$$\begin{matrix}{where} & {{F_{m} = {k^{{- 2}m}\quad{\sum\limits_{n = m}^{\infty}{A_{nm}\left( {k\quad\Delta\quad z} \right)}^{2n}}}},} & {{A_{nm} = {\left( {- 1} \right)^{n - 1}\frac{2^{2n} - 2}{\left( {2n} \right)!}{B_{2n}\begin{pmatrix}n \\m\end{pmatrix}}}},}\end{matrix}$ B_(n) is the n^(th) Bernoulli number, and κ²=k_(x) ²+k_(y)².
 23. An apparatus as claimed in claim 21 and adapted to perform eachsummation over m over the range m=0 to m=1.
 24. An apparatus as claimedin any of claims 17 to 23 and adapted to apply a spatially compactoperator.
 25. An apparatus for determining a vertical component ofparticle velocity from seismic data acquired at first and secondvertically-separated streamers disposed within a water column, the firststreamer comprising M (M>1) pressure sensors disposed along the lengthof the first streamer and the second streamer comprising N (N>1)pressure sensors disposed along the length of the second streamer, theapparatus comprising means for determining the vertical component ofparticle velocity by applying to the pressure acquired at sensors on atleast one streamer an operator having a length (in the spatial domain)that is greater than one but that is less than the number of pressuresensors on the respective streamer.
 26. An apparatus as claimed in anyof claims 17 to 25 and further adapted to process the seismic data usingthe vertical component of the particle velocity.
 27. An apparatus asclaimed in claim 26 and adapted to process the seismic data so as todetermine at least one of an up-going component and a down-goingcomponent of the seismic data.
 28. An apparatus as claimed in claim 26and adapted to process the seismic data using the vertical component ofbase particle velocity thereby to attenuate effects in the processeddata of seismic energy reflected and/or scattered at the free-surface ofthe water column.
 29. An apparatus as claimed in claim 26 and adapted toprocess the seismic data using the vertical component of the particlevelocity thereby to obtain information about the signature of the sourceof seismic energy.
 30. An apparatus as claimed in any of claims 17 to29, and comprising a programmable data processor.
 31. A storage mediumcontaining a program for the data processor of an apparatus as defied inclaim
 30. 32. A storage medium containing a program for controlling adata processor to perform a method as defined in any of claims 1 to 16.